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INTRODUCTION 


In  this  report,  a  long  free-free  slender  beam  is  used  as  a  model  for  a 
flexible  missile  or  rocket.  The  beam  behaves  as  a  Bernoul 1 i-Euler  column,  and 
in  this  case  is  assumed  to  be  rotating  about  its  longitudinal  axis  and  subject 
to  an  end  thrust  (Figure  1).  Of  prime  interest  is  the  effect  of  the  rotation  on 
the  lateral  stability  of  the  beam.  The  motion  is  assumed  to  be  planar. 

Different  phases  of  the  problem  have  been  investigated  in  the  past.  A  sum¬ 
mary  of  the  previous  work  is  given  in  Reference  1.  Silverberg  (ref  2)  was  the 
first  to  include  thrust  on  the  flying  column.  The  differential  equation  for  a 
free-flying  beam  was  given  earlier  as  shown  in  Reference  3.  Beal  (ref  d)  and 
Feodos'ev  (ref  5)  obtained  results  with  pulsating  thrust.  In  1972,  Matsumoto 
and  Mote  (ref  6)  treated  a  similar  problem  with  directional  thrust.  In  this 
case,  however,  feedback  control  was  included  and  a  time  delay  was  applied  to  the 
control.  The  next  contribution  to  understanding  the  problem  was  given  by  Peters 
and  Wu  (ref  1).  They  concentrated  on  mode  shape  solutions  at  zero  frequency  for 
different  thrusts.  A  comorehensi ve  description  is  also  given  in  Reference  1  for 
the  eigenvalues  and  mode  shape  near  zero  thrust  and  with  a  thrust  direction 
dose  to  that  of  a  follower  force.  Recently,  Park  and  Mote  (ref  7)  included  a 
concentrated  mass  and  feedback  control.  The  feedback  control  included  was 
allowed  to  be  from  different  points  along  the  beam. 

As  stated  above,  this  report  assesses  the  effect  of  rotation  on  the  stabil¬ 
ity  of  a  free-free  beam.  The  following  section  is  a  description  of  the  problem. 
Then  the  variational  statement  used  for  the  solution  is  described.  Next  we  show 
how  the  variational  statement  is  used  with  finite  elements  to  solve  the  problem, 
and  lastly,  we  discuss  the  results  of  our  investigation. 

References  are  listed  at  the  end  of  this  report. 


1 


PROBLEM  STATEMENT 


The  geometry  of  the  problem  is  shown  in  Figure  1.  The  beam  has  a  constant 

cross-section  of  area  A,  density  p,  Young's  modulus  E,  and  moment  of  inertia  I. 
It  shows  a  free-flying  column  subject  to  axial  thrust  with  directional  control 
and  rotating  about  its  axis.  The  differential  equation  for  the  beam  is  given  by 

Elu'iv  +  P(|  u '  )  '  +  pAu  +  pAQ2u  =  0  (1) 

The  first  three  terms  represent  the  column  as  treated  in  Reference  1.  The  last 
term  on  the  left-hand  side  shows  the  effect  of  the  rotation.  The  boundary  con¬ 
ditions  are  given  by 

u"  ( 0 )  =  0  ,  u"(  n  =  0  .  u  = 

dr* 

u“  ‘  (0 )  =  0  ,  Elu" ' ( i )  -  KgPu 'll)  -  0  (2) 

In  dimensionless  form  with 

u  =  u/  {  ,  x  -  x/ *  r  -  t /T 


pAf4 

'EI~ 


(V  =  QT 


(3  ) 


and  writing 


U  (  X  .  T  ) 

the  differential  equation  then  becomes 


u ( x le*1 


u'”'  +  Q(xu')'  1  1  u)iU  -  0 


with  the  boundary  conditions 

u" ( 0 )  =  0 
u'"  (0)  =  0 
u" ( 1 )  =  0 


u"'  (1)  -  KfjQfu’  (1)1  i  0 
Rewriting  Eq.  (5)  as  (and  dropping  hats) 

u""  +  Q(xu')'  i  (X2+oj2)u  --  0 


;  1  ) 


h  ) 


(  6  l 


(  T  ) 


It  appears  that  the  addition  of  rotation  simply  shifts  the  frequency  of  vibra- 


tion  of  the  system.  The  boundary  conditions,  Eq.  (6),  become 

u" ( 0 )  =  0 
u‘"  (0)  *  0 
u" ( 1 )  =  0 

u'"(l)  -  KqQu  '  ( 1 )  =  0  (8) 

The  spacial  variables  are  made  dimensionless  by  dividing  through  by  the  beam's 
length  i  and  time  is  made  dimensionless  by  dividing  through  by  a  constant  T  = 
(pAf4,EI)^  which  has  the  units  of  time. 

The  parameter  X  is  a  complex  number  m  general 

X  =  Xp  +  i X  j 

where  both  Xp  and  Xj  are  real  numbers. 

VARIATIONAL  STATEMENT 

To  find  the  form  of  the  variational  statement,  the  differential  eauat-on  - s 
multiplied  by  an  arbitrary  variation  of  the  adjoint  field  variable,  Sv(x),  and 
integrated  over  the  beam  length.  Integrat lon-by-parts  indicates  the  form  of  the 
variational  statement  and  the  natural  boundary  conditions.  The  variational 
statement  is  given  by 

6  J  =  0  I  9  ) 

where 

,1 

J  =  /  [u"v"  -  Qxu ’ v '  +  (Xi+w*  )uv]dx  +  Q( 1  +  Ky  )u ' ( 1 ) v( 1 )  (101 

Derforming  the  variation  of  J  with  respect  to  u  and  v,  one  can  arrive  at  the 
original  boundary  value  problem  as  well  as  the  adjoint.  Eouation  (10)  is  the 
basis  for  a  finite  element  solution  to  the  described  problem. 
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FINITE  ELEMENT  AND  NUMERICAL  FORMULATION 

The  procedure  begins  by  taking  the  variation  of  Eq.  (10)  and  allowing  the 
variations  in  the  problem  variable,  6u(x),  to  be  zero,  i.e.,  varying  adjoint 
variable  v(x)  only  for  now, 


/  [u"6v"  -  Qxu  '  5v  '  +  A2u<5v]dx  -  Q ( 1  +  K0  )  u  '  ( 1 )  6v  ( 1  )  =  0 


<11 


where  A2  =  A2  +  u2 .  To  discretize,  the  beam  is  divided  into  L  elements,  letting 


i-1 

Z  =  L{x - }  i  =  1,2,3 - L  112! 

L 

be  the  running  coordinate  in  each  element.  Substituting  Eq.  (12)  into  Eq.  (Hi 


L 


3  J1  [L3i  M  )"5v< i  )"  -  Q { ^  +  ( i-1 )  (u< 1 ) '6v< 1 ) ' 
L  0 
i  =  l 


-  -  u  < 1 ) 6v  <  ^  ] ds 
L 


-  Q ( 1  +  K0 ) u ( L  > ' ( 1 ) <  L )  (1)  =  0  (13) 

In  order  that  the  displacements  and  their  derivatives  within  an  element  oe 
exoressed  in  terms  of  their  nodal  values,  the  coordinate  vectors  are  introduced. 


Gt1)1  =  jujd) 

U2  <  1  ) 

u3(1) 

u4  M  M 

y(  i  )  T  =  { v 1 < 1 ) 

v2  <  1  ) 

V3M) 

v  4  (  1  >  f 

;  14 

Ujt1),  U2<'')  represent  the  displacement  and  slope  at  the  left  end  of  the  ith 
element,  and  U 3 ( 1 1  and  represent  deflection  and  slope  at  the  right  end.  A 

similar  mterpretat ion  is  applied  to  the  adjoint  coordinate  vector  V ( 1 ) .  The 
transform  is  indicated  by  T. 

Hermitian  polynomials  are  used  to  relate  the  displacements  within  an  ele¬ 
ment  to  its  nodal  values,  hence,  the  following  shape  function  is  assumed: 

aTm  =  jl  -  3H  ♦  2*’  ,  {  -  2£  ’  ♦  V  ,  3£  2  -  2?’  ,  -£2  ♦  Z 1  i  (IS) 

so  that 
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'  16) 


u(  i  )  (C  )  =  aT(? )U< 1  ) 

v(i)(C)  =  aT(^)v('i) 

Substituting  Eq.  (16)  into  Eq.  (13) 

L 

^  D  (  ■>  )  T  {  L  3C  -  Q  [D+  ( i  -  1  )  8  ]  +  --  A }  <5V  ( ’  )  -  Q  [  1  +  Kf)  ]U  t  L  >  TE6V  ( L  >  =  0  (17) 

i  =  1 


,1  -  . 


A  =  /  aaTd£  ,  B  =  /  a'a'Td£  ,  C  =  /  a"a"Td£ 


0  =  /  £a '  a '  Td^  ,  E  =  a* ( L ) aT ( L  ) 
0 


Rewr i t i ng  Eq .  (17), 

L 

^  u(  1  )T{A2P(  1  )  4-  S<  i  )  |6V(  i  )  =  0  (19) 

i  =  1 


where 


p(  i  )  = 

A/L 

i  =  1,2, 

.  .  .  ,  L 

S  ( i'  )  = 

L  3C  - 

Q[D  + 

(i-l)B] 

i  =  1,2, 

... , L- 

1 

S<L)  = 

L'C  - 

Q[D  + 

(L-l)B]  -  Q( 1  +  Kg  )E 

using  certain  continuity  conditions 

between  the 

e  1  ement 

nodal 

va 1 ues 

(i 

U1 

)  U 

=  u3 

-1) 

(i) 

Vl 

(i-l) 

=  v3 

(  i 
u2 

)  (i 

=  U4 

-1) 

(i) 

V2 

( i-l ) 

=  v4 

One  can  write 

-T  ,  U) 

UT  =  {Uj 

(1) 

U2 

(1) 

u3 

(1) 

U4 

(2) 

u3 

(2) 

U4  .  .  .  , 

(L) 

•  u3 

(L)  . 
U4  ) 

-t  ,  (D 
vT  =  {v: 

(1) 

v2 

(1) 

v3 

(I) 

v4 

(2) 

v3 

(2) 

v4  •  •  .  • 

(  L  ) 

•  v3 

(L). 

v4  ) 

(20) 


(  21 ) 


(22) 
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Finally,  [P]  and  [S]  are  NxN  matrices  with  N  =  2L  +  2 .  Since  5v  is  arbitrary,  tne 
eigenvalue  problem  reduces  to 

Ut|A2[P]  +  [S] }  =  0  (23) 

which  is  solved  for  the  eigenvalues. 

CONCLUSIONS  AND  DISCUSSION 

In  this  report,  we  have  included  rotation  about  the  longitudinal  axis  in 
the  dynamic  stability  study  of  a  free-flying  missile  subjected  to  axial  thrusts. 
It  is  assumed  that  the  motions  of  bending  and  the  thrust  are  in  the  same  plane. 
In  the  differential  eguation,  the  only  difference  resulting  from  the  introduc¬ 
tion  of  rotation  is  a  change  in  the  frequency  parameter  A2  to 

A2  =  A2  +  w2  (24) 

where  w  is  the  rotation.  Consequently,  all  the  stability  curves  obtained  pre¬ 
viously  (ref  1)  can  be  used  with  some  simple  modifications.  It  should  be  noted 
that  in  Reference  1,  we  have  written  (with  o  =  0) 

A  =  A=:Ap+iAj  (25) 

and  the  stability  character  of  the  oroblem  is  indicated  by:  (1)  stable  vibra¬ 
tions  =  At  *  0,  A p  =  0;  (2)  unstable  by  buckling  (divergence)  =  Ar  *  0,  Aj  =  0; 
(3)  unstable  by  flutter  =  Ar  *  0,  Aj  *  0;  and  (4)  marginally  stable  =  Aj  =  Ar  = 

0 . 

For  the  present  case,  the  stability  behavior  is  indicated  as  above,  but 
with  Aj  and  Ar  replacing  Aj  and  Ar  in  the  previous  stability  curves 

A=Ar+iAj  (26) 

and 

A 2  =  ( Ar +  i  A t  )  2  =  A 2  +  w-  =  ( Ar ♦  i  A j  )  2  +  w'  (27) 
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or 

A2  =  (Ap+iA;j)2  =  A2  -  w2  =  (Ar+iAt)2  -  w2  ( 28) 

from  Eq.  (28),  when  Ar  =  0,  A2  =  -Aj2  -  u2  ,  hence  Ar  =  0  and  Aj2  =  Aj2+ur‘. 

Thus,  originally  stable  vibrations  will  remain  stable  with  higher  vibration  fre¬ 
quency.  On  the  other  hand,  when  Aj  =  0,  A2  =  Ar2  -  w2 ,  hence  A2  =  Ap2  -  tv2  . 
Thus,  originally  divergent  motions  will  become  stable  vibrations  when  Ap2  <  . 

In  the  case  of  marginal  stability  A  =  0  will  certainly  be  stabilized  since  aj2  = 

U)'  . 

In  the  case  of  flutter  instability,  Eq.  (28)  states  that  A  is  comp'lex  ;A*  * 
0,  Ap  *  0)  if  and  only  if  A  is  complex  ( A j  *  0,  Ar  t  0).  Therefore,  the  f’utter 
instability  is  not  affected  by  the  introduction  of  the  rotation,  which  is  an 
-nteresting  observation. 

Several  demonstrative  stability  curves  with  A2  (and  A2)  versus  Q/n2  are 
shown  in  Eigures  2  through  5.  Only  the  lowest  eigenvalue’s  branches  are  shown, 
Since  they  are  the  ones  which  dictate  the  stability  behavior.  Figure  2  shows 
the  two  lowest  stao1e  vibration  modes  and  two  rigid  body  modes  on  the  A2  -  0 
axis.  T h i s  is  the  case  o*  a  free-flying  missile  with  a  follower  thrust  (Ky  =  C) 
and  with  a  dimensionless  ■'otation  of  w2  =  500.  The  two  fluxural  modes  coalesce 
at  load  Q/n2  =  11.18  beyond  which  flutter  instability  begins.  The  rigid  body 
modes  without  rotation  indicate  margma’’  stab’lity.  Due  to  the  rotation  w,  the 
ax ’S  is  shifted  from  A2  =  0  to  a  2  =  0,  therefore,  these  previously  rigid  body 
modes  are  now  stable  modes  of  vibrations.  The  thrust  that  'S  controlled  with  a 
small  negative  tangency  (Ky  =  -0.05)  is  shown  in  Figure  3.  It  is  noted  in  this 
figure  that  the  divergence  instability  without  rotation  is  stabilized  by  u>2  = 
500.  However,  the  new  critical  load  is  lowered  from  Q/n2  =  11.18  to  5.30,  not 
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because  of  u*  ,  but  due  to  the  negative  control  parameter  Kg.  Figure  4  shows 
case  of  Ky  =  -l  or  that  the  thrust  has  a  fixed  direction  of  the  inertia  ax>s. 
It  is  clear  that  the  divergence  instability  of  the  lowest  branch  is  stabilize! 
so  that  the  critical  load  has  been  raised  from  zero  to  Qqr  =  1.50  it'.  Fm.y  ' 
the  case  for  a  small  positive  tangency  control  parameter  (Kg  =  0.05)  -is  shown 
figure  5.  In  this  figure,  the  original  divergence  instability  at  Q ./ rr ■ '  =  3.00 
stabilized  by  u' .  However,  the  original  critical  load  of  flutter  instability 
Q  n:  =  9.90  is  not  changed  by  the  rotation.  Hence,  the  critical  load  in  tfrs 
case  ’s  raised  from  3.00  to  9.90  due  to  the  rotation  of  w:  =  500. 
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